

clear all
set more off
cap log close

***************************************************************************************************
* Program: analysis.do
* Purpose: Produce tables/figures
***************************************************************************************************

************************************************************************
** set global macros
************************************************************************
global root_dir "~/Dropbox/Prison_Covid/Visitation misconduct paper/Science Advances R&R/Replication"
global data_raw "${root_dir}/Data/Raw"
global gas_price_dir "${data_raw}/gas_price"
global compstat_dir "${data_raw}/cdcr_compstat"
global macroecon_data "${data_raw}/macro_economic_data"
global boundaries_dir "${data_raw}/boundaries"
global fire_dir "${data_raw}/california_historic_fire_perimeters_data"
global weather_dir "${data_raw}/weather"
global processed_dir "${root_dir}/Data/Processed"
global output "${root_dir}/output"

*****************************************************************************
* Figure 1: gas price and prison incident (by year)
*******************************************************************************

use "${processed_dir}/CA_incidents_visit_reg_data.dta", clear 

collapse (mean) monthly_gas_price incidents_per100, by(year month)

gen date = mdy(month, 1, year)
format date %td
drop if month > 3 & year == 2020


sort year month 
tw (connected monthly_gas_price date if date < mdy(4,1,2020), msize(vsmall) yaxis(1) ytitle("Gas Price", axis(1))) ///
	(connected incidents_per100 date if date < mdy(4,1,2020), msize(vsmall) yaxis(2) ytitle("Incidents Per 100", axis(2)) ///
	legend(order(1 "Gas Price" 2 "Incidents Per 100")) ///
	xlabel(`=mdy(1,1,2017)'(180)`=mdy(4,1,2020)', format(%tdMonCCYY) labsize(small)) ///
	xtitle("Date"))
graph export "${output}/figure1.jpg", replace width(4000)


****************************************************************
* Table 1: summary stat tables 
****************************************************************


use "${processed_dir}/CA_incidents_visit_reg_data.dta", clear 

egen tag_prison = tag(institution)

* Calculate mean distance
quietly summarize wavg_dist_zcta_to_prison if ~missing(norm_visits),
scalar mean_dist = r(mean)

* Create indicators for below and above mean distance
gen below_mean_dist = wavg_dist_zcta_to_prison < mean_dist
gen above_mean_dist = wavg_dist_zcta_to_prison >= mean_dist


label variable tot_visits_7to15 "Total Visits"
label variable norm_visits_7to15 "Normalized Visits"

****************************************************************

* Clear previous estimates
eststo clear
estimates clear

* Whole sample
eststo whole: quietly estpost summarize ///
	tot_visits_7to15 norm_visits_7to15 wavg_dist_zcta_to_prison ///
	incidents_per100 uof_per100 battery assault drug ///
	appeals_per100 disciplinaries_per100 placement_csw if ~missing(norm_visits), detail

* Below mean distance
eststo below: quietly estpost summarize ///
	tot_visits_7to15 norm_visits_7to15 wavg_dist_zcta_to_prison ///
	incidents_per100 uof_per100 battery assault drug ///
	appeals_per100 disciplinaries_per100 placement_csw if below_mean_dist & ~missing(norm_visits), detail

* Above mean distance
eststo above: quietly estpost summarize ///
	tot_visits_7to15 norm_visits_7to15 wavg_dist_zcta_to_prison ///
	incidents_per100 uof_per100 battery assault drug ///
	appeals_per100 disciplinaries_per100 placement_csw if above_mean_dist & ~missing(norm_visits), detail

* Calculate diff	
eststo diff: estpost ttest ///
	tot_visits_7to15 norm_visits_7to15 wavg_dist_zcta_to_prison ///
	incidents_per100 uof_per100 battery assault drug ///
	appeals_per100 disciplinaries_per100 placement_csw ///
	if (below_mean_dist | above_mean_dist) & ~missing(norm_visits), ///
	by(above_mean_dist)
	
* Export summary stats
esttab whole below above diff using ///
	"${output}/table1.tex", ///
	label tex replace mtitles("Full Sample" "Below Mean Distance" "Above Mean Distance" "Diff") ///
	cells("mean(pattern(1 1 1 0) fmt(3)) sd(pattern(1 1 1 0)) b(star pattern(0 0 0 1) fmt(3)) t(pattern(0 0 0 1) par fmt(3))")
	
******************************************************************************************
* MAIN ANALYSIS
* Monthly Regression anaylsis 
******************************************************************************************

use "${processed_dir}/CA_incidents_visit_reg_data.dta", clear 

global all_ctrls "tmean tmax ppt ca_unemployment_rate log_sp500 log_overtime_hours log_sick_leave_hours wildfire"
global partial_ctrls "tmean tmax ppt ca_unemployment_rate log_sp500 wildfire"

drop if mid_date >= mdy(4, 1, 2020)

********************************
* Asymmetric Residuals
********************************
reg log_monthly_gas_price i.month tmean tmax ppt, robust
predict resid, resid
* Asymmetric effects 
gen resid_pos = resid * (resid > 0)
gen resid_neg = -1 * resid * (resid < 0)
//gen resid_neg_abs = -1 * resid * (resid < 0)

label variable resid_pos "Gas price residual X 1 X (residual $>$ 0)"
label variable resid_neg "Gas price residual X -1 X(residual $<$ 0)"
		
********************************
* Table 2 
********************************
estimates clear 
eststo clear 
* OLS 
eststo ols: reghdfe incidents_per100 log_norm_visits_7to15 $partial_ctrls if ~missing(log_norm_visits_7to15), cluster(institution) absorb(institution month) 
estadd local prisonfe "Yes"
estadd local monthfe "Yes"

* 2sls
eststo tsls: ivreghdfe incidents_per100 $partial_ctrls ///
	(log_norm_visits_7to15=log_monthly_gas_price) if ~missing(log_norm_visits_7to15), ///
	cluster(institution) absorb(institution month) first
estadd local prisonfe "Yes"
estadd local monthfe "Yes"
* Add Kleibergen–Paap, Cragg–Donald, and Hansen only where available
capture local tmp1 = round(`=e(idstat)', .01)
capture estadd local elmstat "`tmp1'"
capture local tmp2 = round(`=e(widstat)', .01)
capture estadd local ewidstat "`tmp2'"


* first stage (Base)
eststo fs1: reghdfe log_norm_visits_7to15 log_monthly_gas_price $partial_ctrls if ~missing(log_norm_visits_7to15), ///
	cluster(institution) absorb(institution month) 
estadd local prisonfe "Yes"  
estadd local monthfe "Yes"

* First stage: visitation on unexpected gas shocks (Asymmetric)
eststo fs2: reghdfe log_norm_visits_7to15 resid_pos resid_neg ca_unemployment_rate log_sp500 wildfire if ~missing(log_norm_visits_7to15), cluster(institution) absorb(institution month)
estadd local prisonfe "Yes"
estadd local monthfe "Yes"

* reduce form
eststo rf1: reghdfe incidents_per100 log_monthly_gas_price $partial_ctrls if ~missing(log_norm_visits_7to15), ///
	cluster(institution) absorb(institution month) 
estadd local prisonfe "Yes"  
estadd local monthfe "Yes"

* Reduced form (with controls)
eststo rf2: reghdfe incidents_per100 resid_pos resid_neg ca_unemployment_rate log_sp500 wildfire if ~missing(log_norm_visits_7to15), absorb(institution month) cluster(institution)
estadd local prisonfe "Yes"
estadd local monthfe "Yes"

esttab ols tsls fs1 fs2 rf1 rf2 using ///
"${output}/table2.tex", ///
replace label booktabs nonotes collabels(none) ///
b(4) se(4) star(* 0.10 ** 0.05 *** 0.01) ///
mtitles("OLS" "2SLS" ///
        "First Stage" "First Stage" "Reduced Form" "Reduced Form") ///
scalars("prisonfe Prison FE" ///
		"monthfe Month FE" ///
        "elmstat Kleibergen-Paap rk LM statistic (Underidentification test)" ///
        "ewidstat Cragg-Donald Wald F statistic (Weak identification test)") ///
mgroups("Incidents per 100" "Log(Normalized Visits)" "Incidents per 100", ///
        pattern("1 0 1 0 1 0") ///
        prefix(\multicolumn{@span}{c}{) suffix(}) span ///
        erepeat(\cmidrule(lr){@span}))
		
*********************************************
* Main Table 3
*********************************************

eststo clear
estimates clear 

local vars "incidents_per100 uof_per100 battery assault drug appeals_per100 disciplinaries_per100 placement_csw"
foreach var of local vars{
	eststo: ivreghdfe `var' $partial_ctrls (log_norm_visits_7to15=log_monthly_gas_price) if ~missing(log_norm_visits), cluster(institution) absorb(institution month) 
	estadd local prisonfe "Yes" 
	estadd local monthfe "Yes"
	local tmp1 = round(`=e(idstat)',.01)
	estadd local elmstat "`tmp1'"	
	local tmp2 = round(`=e(widstat)',.01)
	estadd local ewidstat "`tmp2'"
// 	local tmp3 = round(`=e(j)',.01)
// 	estadd local ehansen "`tmp3'"
	
}

esttab using "${output}/table3.tex", ///
 replace se b(4) se(4) star(* 0.10 ** 0.05 *** 0.01)   ///
 label booktabs nonotes collabels(none) ///
	scalars("prisonfe Prison FE" "monthfe Month FE" "elmstat Kleibergen-Paap rk LM statistic (Underidentification test)" ///
	"ewidstat Cragg-Donald Wald F statistic (Weak identification test)") 

 
 
********************************
* Table S4: Robustness - First Stage
********************************
	
estimates clear 
eststo clear
local vars "log_norm_visits_7to15" 
foreach var of local vars{
	
	eststo: reghdfe `var' log_monthly_gas_price $partial_ctrls, cluster(institution) absorb(institution month)
	estadd local prisonfe "Yes"
	estadd local monthfe "Yes"
	
	eststo: reghdfe `var' log_monthly_gas_price $partial_ctrls if year < 2020, cluster(institution) absorb(institution month)
	estadd local prisonfe "Yes"
	estadd local monthfe "Yes"
	
	eststo: reghdfe `var' log_monthly_gas_price $all_ctrls , cluster(institution) absorb(institution month)
	estadd local prisonfe "Yes"
	estadd local monthfe "Yes"

	eststo: reghdfe `var' log_monthly_gas_price $partial_ctrls if ~inlist(month, 11, 12, 1) , cluster(institution) absorb(institution month)
	estadd local prisonfe "Yes"
	estadd local monthfe "Yes"

}

local vars "log_norm_visits log_norm_visits_7to15_1hr" 
foreach var of local vars{
	
	eststo: reghdfe `var' log_monthly_gas_price $partial_ctrls , cluster(institution) absorb(institution month)
	estadd local prisonfe "Yes"
	estadd local monthfe "Yes"

}

eststo: reghdfe log_norm_visits_new log_monthly_gas_price $partial_ctrls if (year >= 2018) | (year == 2017 & month >= 8), cluster(institution) absorb(institution month)
estadd local prisonfe "Yes"  
estadd local monthfe "Yes"

eststo: reghdfe log_norm_visits_repeat log_monthly_gas_price $partial_ctrls if (year >= 2018) | (year == 2017 & month >= 8), cluster(institution) absorb(institution month)
estadd local prisonfe "Yes"  
estadd local monthfe "Yes"

local vars "log_norm_visits_night" 
foreach var of local vars{
	
	eststo: reghdfe `var' log_monthly_gas_price $partial_ctrls , cluster(institution) absorb(institution month)
	estadd local prisonfe "Yes"
	estadd local monthfe "Yes"

}

esttab using ///
"${output}/tables4.tex", ///
replace label booktabs nonotes collabels(none) ///
b(4) se(4) star(* 0.10 ** 0.05 *** 0.01) ///
mtitles("Baseline" "2017-2019"  "Add Prison Strain Ctrls (2017-2019)" ///
		"Exclude Winter Months" "Weekend (All Hours)" ///
		"Weekend Visits $>=$1 Hour" ///
		"New Visits" ///
		"Repeat Visits" ///
		"Placebo: 0-4am Presence") ///
scalars("prisonfe Prison FE" ///
		"monthfe Month FE") ///
mgroups("Log(Normalized Visits)", ///
        pattern("1 0 0 0 0 0 0") ///
        prefix(\multicolumn{@span}{c}{) suffix(}) span ///
        erepeat(\cmidrule(lr){@span}))
		
		
********************************
* Table S5: Robustness - 2SLS
********************************

set more off
estimates clear 
eststo clear
local vars "log_norm_visits_7to15" 
foreach var of local vars{
	
	eststo: ivreghdfe incidents_per100 $partial_ctrls ///
		(`var'=log_monthly_gas_price) if ~missing(`var'), ///
		cluster(institution) absorb(institution month) first
	estadd local prisonfe "Yes"
	estadd local monthfe "Yes"
	local tmp1 = round(`=e(idstat)',.01)
	estadd local elmstat "`tmp1'"	
	local tmp2 = round(`=e(widstat)',.01)
	estadd local ewidstat "`tmp2'"
	
	eststo: ivreghdfe incidents_per100 $partial_ctrls ///
		(`var'=log_monthly_gas_price) if ~missing(`var') & year < 2020, ///
		cluster(institution) absorb(institution month) first
	estadd local prisonfe "Yes"
	estadd local monthfe "Yes"
	local tmp1 = round(`=e(idstat)',.01)
	estadd local elmstat "`tmp1'"	
	local tmp2 = round(`=e(widstat)',.01)
	estadd local ewidstat "`tmp2'"
	
	eststo: ivreghdfe incidents_per100 $all_ctrls ///
		(`var'=log_monthly_gas_price) if ~missing(`var'), ///
		cluster(institution) absorb(institution month) first
	estadd local prisonfe "Yes"
	estadd local monthfe "Yes"
	local tmp1 = round(`=e(idstat)',.01)
	estadd local elmstat "`tmp1'"	
	local tmp2 = round(`=e(widstat)',.01)
	estadd local ewidstat "`tmp2'"
	
	eststo: ivreghdfe incidents_per100 $partial_ctrls ///
		(`var'=log_monthly_gas_price) if ~missing(`var') & ~inlist(month, 11, 12, 1), ///
		cluster(institution) absorb(institution month) first
	estadd local prisonfe "Yes"
	estadd local monthfe "Yes"
	local tmp1 = round(`=e(idstat)',.01)
	estadd local elmstat "`tmp1'"	
	local tmp2 = round(`=e(widstat)',.01)
	estadd local ewidstat "`tmp2'"

}

label variable log_norm_visits_7to15_1hr "Log Normalized Visits ($>=$1 Hour)"

local vars "log_norm_visits log_norm_visits_7to15_1hr" 
foreach var of local vars{

	
	eststo: ivreghdfe incidents_per100 $partial_ctrls ///
		(`var'=log_monthly_gas_price) if ~missing(`var') , ///
		cluster(institution) absorb(institution month) first
	estadd local prisonfe "Yes"
	estadd local monthfe "Yes"
	local tmp1 = round(`=e(idstat)',.01)
	estadd local elmstat "`tmp1'"	
	local tmp2 = round(`=e(widstat)',.01)
	estadd local ewidstat "`tmp2'"


}

label variable log_norm_visits_new "Log Normalized Visits (New)"
label variable log_norm_visits_repeat "Log Normalized Visits (Repeat)"

local vars "log_norm_visits_new log_norm_visits_repeat" 
foreach var of local vars{

	
	eststo: ivreghdfe incidents_per100 $partial_ctrls ///
		(`var'=log_monthly_gas_price) if ~missing(`var') & ((year >= 2018) | (year == 2017 & month >= 8)), ///
		cluster(institution) absorb(institution month) first
	estadd local prisonfe "Yes"
	estadd local monthfe "Yes"
	local tmp1 = round(`=e(idstat)',.01)
	estadd local elmstat "`tmp1'"	
	local tmp2 = round(`=e(widstat)',.01)
	estadd local ewidstat "`tmp2'"


}

label variable log_norm_visits_night "Log Normalized Visits (0-4am Presence)"
local vars "log_norm_visits_night" 
foreach var of local vars{

	
	eststo: ivreghdfe incidents_per100 $partial_ctrls ///
		(`var'=log_monthly_gas_price) if ~missing(`var'), ///
		cluster(institution) absorb(institution month) first
	estadd local prisonfe "Yes"
	estadd local monthfe "Yes"
	local tmp1 = round(`=e(idstat)',.01)
	estadd local elmstat "`tmp1'"	
	local tmp2 = round(`=e(widstat)',.01)
	estadd local ewidstat "`tmp2'"


}


esttab using ///
"${output}/tables5.tex", ///
replace label booktabs nonotes collabels(none) ///
b(4) se(4) star(* 0.10 ** 0.05 *** 0.01) ///
mtitles("Baseline"  "2017-2019"  "Add Prison Strain Ctrls (2017-2019)" ///
		"Exclude Winter Months" "Weekend (All Hours)" ///
		"Weekend Visits $>=$1 Hour" ///
		"New Visits" ///
		"Repeat Visits" ///
		"Placebo: 0-4am Presence") ///
scalars("prisonfe Prison FE" ///
		"monthfe Month FE" ///
        "elmstat Kleibergen-Paap rk LM statistic (Underidentification test)" ///
        "ewidstat Cragg-Donald Wald F statistic (Weak identification test)") ///
mgroups("Incidents per 100 Prisoners", ///
        pattern("1 0 0 0 0 0 0 0 0") ///
        prefix(\multicolumn{@span}{c}{) suffix(}) span ///
        erepeat(\cmidrule(lr){@span}))
		
		
	
*********************************************
* Table S6: Using both price and log_ca_priceXdist as instruments
*********************************************
estimates clear 
eststo clear 
* OLS 
eststo: reghdfe incidents_per100 log_norm_visits_7to15 $partial_ctrls if ~missing(log_norm_visits), cluster(institution) absorb(institution month) 
estadd local prisonfe "Yes" 
estadd local monthfe "Yes"
* 2sls: two instruments
eststo: ivreghdfe incidents_per100 $partial_ctrls ///
	(log_norm_visits_7to15=log_monthly_gas_price log_ca_priceXdist) if ~missing(log_norm_visits), ///
	cluster(institution) absorb(institution month) first
estadd local prisonfe "Yes"  
estadd local monthfe "Yes"
local tmp1 = round(`=e(idstat)',.01)
estadd local elmstat "`tmp1'"	
local tmp2 = round(`=e(widstat)',.01)
estadd local ewidstat "`tmp2'"
local tmp3 = round(`=e(j)',.01) //Hansen 
estadd local ehansen "`tmp3'" 
local tmp4 = round(`=e(jp)',.01) // Hansen -p 
estadd local phansen "`tmp4'"
	
* first stage
eststo: reghdfe log_norm_visits_7to15 log_monthly_gas_price log_ca_priceXdist $partial_ctrls if ~missing(log_norm_visits), ///
	cluster(institution) absorb(institution month) 
estadd local prisonfe "Yes"  
estadd local monthfe "Yes"	
* reduce form
eststo: reghdfe incidents_per100 log_monthly_gas_price log_ca_priceXdist $partial_ctrls if ~missing(log_norm_visits), ///
	cluster(institution) absorb(institution month) 
estadd local prisonfe "Yes"  
estadd local monthfe "Yes"	
* export to table	
* OLD table name: reghdfe_log_incident_gas_price_2sls_two_ivs.tex
 esttab using "${output}/tables6.tex", ///
 replace se b(4) se(4) star(* 0.10 ** 0.05 *** 0.01)   ///
 label booktab nonotes collabels(none) mtitles("OLS" "2SLS" "First Stage" "Reduced Form") ///
 scalars("prisonfe Prison FE" "monthfe Month FE" "elmstat Kleibergen-Paap rk LM statistic (Underidentification test)" ///
	"ewidstat Cragg-Donald Wald F statistic (Weak identification test)" ///
	"ehansen Hansen J statistic (Overidentification test)" "phansen p-value for Hansen J statistic") ///
 mgroups("Incidents per 100" "Log(Normalized Visits)" "Incidents per 100", pattern("1 0 0 1") ///
 prefix(\multicolumn{@span}{c}{) suffix(}) span erepeat(\cmidrule(lr){@span}))
 	 
************************************************************************************************
* Table S7: Lead/lag price to understand the duration of effects
************************************************************************************************

 
**********************************************
** Distributed lag model + placebo test ***
**********************************************

estimates clear 
eststo clear 
eststo: reghdfe incidents_per100 $partial_ctrls ///
	log_monthly_gas_price log_gas_price_post1 if ~missing(log_norm_visits), ///
	cluster(institution) absorb(institution month)
estadd local prisonfe "Yes" 
estadd local monthfe "Yes"
eststo: reghdfe incidents_per100 $partial_ctrls ///
	log_monthly_gas_price log_gas_price_post1 log_gas_price_post2 log_gas_price_post3 if ~missing(log_norm_visits), ///
	cluster(institution) absorb(institution month)
estadd local prisonfe "Yes" 
estadd local monthfe "Yes"
eststo: reghdfe incidents_per100 $partial_ctrls ///
	log_monthly_gas_price log_gas_price_pre1 if ~missing(log_norm_visits), ///
	cluster(institution) absorb(institution month)	
estadd local prisonfe "Yes"
estadd local monthfe "Yes" 	
eststo: reghdfe incidents_per100 $partial_ctrls ///
	log_monthly_gas_price log_gas_price_pre1 log_gas_price_pre2 if ~missing(log_norm_visits), ///
	cluster(institution) absorb(institution month)
estadd local prisonfe "Yes"
estadd local monthfe "Yes" 	
eststo: reghdfe incidents_per100 $partial_ctrls ///
	log_monthly_gas_price log_gas_price_pre1 log_gas_price_pre2 log_gas_price_pre3 ///
	log_gas_price_pre4 log_gas_price_pre5 ///
	if ~missing(log_norm_visits), ///
	cluster(institution) absorb(institution month)	
estadd local prisonfe "Yes"
estadd local monthfe "Yes" 
esttab using "${output}/tables7.tex", ///
 replace se b(4) se(4) star(* 0.10 ** 0.05 *** 0.01)   ///
  scalars("prisonfe Prison FE" "monthfe Month FE") ///
 label booktabs nonotes collabels(none) nomtitles ///
 mgroups("Incidents per 100", pattern("1 0 0 0 0") ///
 prefix(\multicolumn{@span}{c}{) suffix(}) span erepeat(\cmidrule(lr){@span}))
 
************************************************************************************************
* Table S8
************************************************************************************************

eststo clear
estimates clear 

local vars "agg_battery agg_battery_nonpri agg_battery_officer assault_on_inmate assault_on_offi_nonpri battery_on_inmate battery_on_offi_nonpri disciplinaries_fighting disciplinaries_cellphone"
foreach var of local vars{
	eststo: ivreghdfe `var' $partial_ctrls (log_norm_visits_7to15=log_monthly_gas_price) if ~missing(log_norm_visits), cluster(institution) absorb(institution month) 
	estadd local prisonfe "Yes" 
	estadd local monthfe "Yes"
	local tmp1 = round(`=e(idstat)',.01)
	estadd local elmstat "`tmp1'"	
	local tmp2 = round(`=e(widstat)',.01)
	estadd local ewidstat "`tmp2'"
// 	local tmp3 = round(`=e(j)',.001)
// 	estadd local ehansen "`tmp3'"
	
}


esttab using "${output}/tables8.tex", ///
 replace se b(4) se(4) star(* 0.10 ** 0.05 *** 0.01)   ///
 label booktabs nonotes collabels(none) ///
	scalars("prisonfe Prison FE" "monthfe Month FE" "elmstat Kleibergen-Paap rk LM statistic (Underidentification test)" ///
	"ewidstat Cragg-Donald Wald F statistic (Weak identification test)") 
	 
	 
********************************
* Appendix Figure S1
********************************
		
cd "${output}"
	
eststo clear 
estimates clear 
forval i = 1/4{
eststo a`i': reghdfe log_norm_visits_7to15 log_monthly_gas_price if dist_quartile == `i', cluster(institution) absorb(institution month) 
eststo b`i': reghdfe log_norm_visits_7to15 log_monthly_gas_price $partial_ctrls if dist_quartile == `i', cluster(institution) absorb(institution month) 
}

coefplot (a*, msize(small) pstyle(p1)) ///
		 (b*, msize(small) pstyle(p2)), ///
 keep(log_monthly_gas_price) asequation swapnames  vertical yline(0) ///
	 coeflabels(, wrap(19)) scheme(s2color) graphregion(color(white)) ///
	 mlabel format(%9.3g) mlabposition(12) mlabgap(*1) ///
	 legend(order(2 "No Controls" 4 "Full Controls") rows(1) size(*0.95)) ///
	 title(, size(medium)) yscale(titlegap(3)) xtitle("") title("Panel (A)") ///
	 eqrename(a1 = "First Quartile" a2 = "Second Quartile" a3 = "Third Quartile" a4 = "Fourth Quartile" ///
	 b1 = "First Quartile" b2 = "Second Quartile" b3 = "Third Quartile" b4 = "Fourth Quartile") name(first_stage, replace)
 
**
eststo clear 
estimates clear 
forval i = 1/4{
eststo a`i': reghdfe incidents_per100 log_monthly_gas_price if ~missing(log_norm_visits) & dist_quartile == `i', cluster(institution) absorb(institution month) 
eststo b`i': reghdfe incidents_per100 log_monthly_gas_price $partial_ctrls if ~missing(log_norm_visits) & dist_quartile == `i', cluster(institution) absorb(institution month) 
}

coefplot (a*, msize(small) pstyle(p1)) ///
		 (b*, msize(small) pstyle(p2)), ///
 keep(log_monthly_gas_price) asequation swapnames  vertical yline(0) ///
	 coeflabels(, wrap(19)) scheme(s2color) graphregion(color(white)) ///
	 mlabel format(%9.3g) mlabposition(12) mlabgap(*1) ///
	 legend(order(2 "No Controls" 4 "Full Controls") rows(1) size(*0.95)) ///
	 title(, size(medium)) yscale(titlegap(3)) xtitle("") title("Panel (B)") ///
	 eqrename(a1 = "First Quartile" a2 = "Second Quartile" a3 = "Third Quartile" a4 = "Fourth Quartile" ///
	 b1 = "First Quartile" b2 = "Second Quartile" b3 = "Third Quartile" b4 = "Fourth Quartile") name(reduced_form, replace)

graph combine first_stage reduced_form, ///
    col(1) ///
    xsize(6) ///
    ysize(9) ///
    imargin(small) ///
    graphregion(color(white))
graph export "${output}/figures1.jpg", ///
    width(4000) replace
	
******************************************************************************************
*  Table S1
******************************************************************************************

use "${processed_dir}/wavg_dist_zcta_to_prison.dta", clear

sort dist_rank

keep prison_name wavg_dist_zcta_to_prison dist_rank

export delimited "${output}/tables1.csv", replace


******************************************************************************************
*  Table S2
******************************************************************************************

use "${processed_dir}/weekly_state_prisons_analysis_data_2017_2020.dta", clear

eststo clear 
estimates clear 
eststo s1: quietly estpost summarize avg_income_1k_vh avg_pct_hisp_vh avg_pct_white_vh avg_pct_black_vh avg_pct_asian_vh avg_pct_college_vh, detail
esttab s1 using ///
	"${output}/tables2.tex", /// 
	label tex replace cells("mean(fmt(3)) sd(fmt(3)) p5(fmt(3)) p25(fmt(3)) p50(fmt(3)) p75(fmt(3)) p95(fmt(3))")
	
	
	
******************************************************************************************
*Table S3
******************************************************************************************

use "${processed_dir}/weekly_state_prisons_analysis_data_2017_2020.dta", clear
 
estimates clear 
eststo clear
local vars "avg_income_1k avg_pct_hisp avg_pct_white avg_pct_black avg_pct_asian avg_pct_college"
foreach var of local vars{
	eststo: reghdfe `var'_vh log_gas_price  public_holiday tmax tmin ppt log_sp500 ca_unemployment_rate, ///
	absorb(Shapefile_ID month) cluster(Shapefile_ID)
	estadd local prisonfe "Yes"
	estadd local monthfe "Yes"
}	
esttab using "${output}/tables3.tex", ///
 replace se b(4) se(4) r2 star(* 0.10 ** 0.05 *** 0.01)   ///
 label booktabs nonotes collabels(none)  ///
 mtitles("Median Household Income (1K)" "\% Hispanic" "\% White" "\% Black" "\% Asian" "\% College") ///
 scalars("prisonfe Prison FE" "monthfe Month FE") ///
  mgroups("Average Demographics for Weekend 7am-3pm Visitors", pattern("1 0 0 0 0 0") ///
 prefix(\multicolumn{@span}{c}{) suffix(}) span erepeat(\cmidrule(lr){@span}))
 
 